!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! Coded by: B. G.
! Date    : 2016-06-28
! Name    : Atmosphere.f95
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! Version:
!    0.0   2016-06-28   BG   Create
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! This routine calculates atmospheric properties of the 
!    US Standard Atmosphere, 1962. Assume an inverse 
!    square gravitational field. this assumption yields 
!    data that agrees with the COESA document within 1%
!    at all altitudes up to 700 kilometers (2296588 ft).
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! Caution:
!   input and output are ISO units, but US units used in 
!    the box.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! Input:
!    H   height, m
! Output:
!    p   pressure, Pa
!    T   temperature, K
!    a   speed of sound, m/s
!    rho density, kg/m**3
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!

program Atmosphere
   !variables declaration
   implicit none
   real*8::H
   real*8::p,T,a,rho
   real*8::G0,Wm0,R0,GMRS
   real*8::G,Z
   real*8::ELH,TMS,ELZ,DMDZ,EM
   real*8::ZLZ
   real*8,dimension(10)::HG
   real*8,dimension(14)::HM
   real*8,dimension(14)::WM
   real*8,dimension(23)::TM
   real*8,dimension(22)::PM
   integer*4::errCode
   integer::i,j,k
   
   data G0,Wm0,R0,GMRS/32.1740484,28.9644,20890855.0,0.018743428/
   data HG/-16404.,0.0,36089.,65617.,104987.,154199.,170604.,200131.,&
         &259186.,291160./
   data HM/295276.,328084.,360892.,393701.,492126.,524934.,557743.,623360.,&
         &754593.,984252.,1312336.,1640420.,1968504.,2296588./
   data WM/28.9644,28.88,28.56,28.07,26.92,26.66,26.4,25.85,24.7,22.66,19.94,&
         &17.94,16.84,16.17/
   data TM/577.17,518.67,389.97,389.97,411.57,487.17,487.17,454.77,325.17,&
         &325.17,379.17,469.17,649.17,1729.17,1999.17,2179.17,2431.17,2791.17,&
         &3295.17,3889.17,4357.17,4663.17,4861.17/
   data PM/3711.0839,2116.2165,472.67563,114.34314,18.128355,2.3162178,&
         &1.2321972,3.8030279E-01,2.1671352E-02,3.4313478E-03,6.2773411E-04,&
         &1.5349091E-04,5.2624212E-05,1.0561806E-05,7.7083076E-06,5.8267151E-06,&
         &3.5159854E-06,1.4520255E-06,3.9290563E-07,8.4030242E-08,2.2835256E-08,&
         &7.1875452E-09/
   
   !calculation process
   !proccess start
   !input height
   write(*,*) 'Please Input Height in meter (0 <= H <= 70 0000 m):'
   read(*,*,iostat=errCode) H
   if ((errCode /= 0) .or. (H < 0) .or. (H > 700000)) then   !input error
      write(*,*) 'Input Error.'
      stop
   end if
   
   H=H*3.28083333   !transform into ft
   
   G=G0*(R0/(R0+H))**2
   if (H <= 295276.0) then
      Z=R0*H/(R0+H)
      do i=2,10
         j=i-1
         if (HG(i) >= Z) exit
      end do
      ELH=(TM(j+1)-TM(j))/(HG(j+1)-HG(j))
      TMS=TM(j)+ELH*(Z-HG(i))
      ELZ=ELH*G/G0
      DMDZ=0.0
      EM=Wm0
      if (ELH == 0.0) then
         p=PM(j)*(TM(j)/TMS)**(GMRS/ELH)
      else
         p=PM(j)*EXP(GMRS*(HG(j)-Z)/TMS)
      end if
   else
      do i=2,14
         j=i+8
         k=i-1
         if (HM(i) >= H) exit
      end do
      ELZ=(TM(j+1)-TM(j))/(HM(k+1)-HM(k))
      TMS=TM(j)+ELZ*(H-HM(k))
      DMDZ=(WM(k+1)-WM(k))/(HM(k+1)-HM(k))
      EM=WM(k)+DMDZ*(H-HM(k))
      ZLZ=Z-TMS/ELZ
      p=PM(j)*EXP(GMRS/ELZ*(R0/(R0+ZLZ))**2*((H-HM(k))*(R0+ZLZ)/((R0+H)*(R0+HM(k)))-&
         &ALOG(TMS*(R0+HM(k))/(TM(j)*(R0+H)))))
   end if
   a=49.022164*SQRT(TMS)
   rho=GMRS*p/(G0*TMS)
   T=EM*TMS/Wm0
   !unit transform
   H=H/3.28083333
   p=p*47.88025898
   a=a*0.3048
   rho=rho*515.378818393
   T=T/1.8
   !output
   write(*,*) 'THE RESULTS ARE:'
   write(*,*) '  H(m)               p(Pa)             a(m/s)            rho(kg/m**3)      T(K)'
   write(*,*) H,p,a,rho,T
   !process compelte
end program Atmosphere
